
rm(list=ls())
setwd("C:/Users/patri/Dropbox/My Projects/Spring 2021/Empirical IR/Data")
d <- read.csv("threatBJPS.csv")

library(MCMCpack)
library(doBy)
set.seed(03112021)

row.names(d) <- d$cntry_yr



d$strivalry_dum2 <- ordered(as.factor(d$strivalry_dum2))
d$hadispute2 <- ordered(as.factor(d$hadispute))
d$bufferstate2 <- ordered(as.factor(d$bufferstate2))


post_us <- MCMCmixfactanal(~ strivalry_dum2+  hadispute2 + std_lntotborder
                            + bufferstate2 +std_meansipri_milgdp3 +std_lntyrs_2,,
                           factors=1, data=d,
                           lambda.constraints = list(strivalry_dum2=list(2,"+")),
                           burnin=20000, mcmc=100000, thin=50,
                           verbose=0, L0=.25, store.lambda=TRUE,
                           store.scores=TRUE, tune=1.2)


save(post_us,file="IRT_Threat_BJPS.Rda")
load("IRT_Threat_BJPS.Rda")



summary(post_us)

# convert object from mcmc to data.frame and transpose
c<-as.data.frame(post_us)
c<-t(c)
c1<-apply(c,1,quantile,probs = c(0.025, 0.05, 0.5, 0.95, 0.975))
mean<-apply(c,1,mean)
sd<-apply(c,1,sd)
c1<-rbind(c1,mean,sd)
c1<-t(c1)

# isolate and extract lambda/psi (beginning and end of variables)
print(c1[11,])
lambda<-c1[1:11,]
n<-nrow(c1)
psi<-c1[(n-2):n,]
lambda<-rbind(lambda,psi)

# isolate country ideal points (phi; in middle)
phi<-c1[12:(n-3),]

# change vars to lower 95+ 90+ median+ and upper 90+ 95; change names to remove ".2"
lambda
colnames(lambda) <- c("low95","low90","median","high90","high95","mean","sd")

# change vars to lower 95+ median+ and upper 90; change names to leave just country name
phi
colnames(phi) <- c("low95","low90","median","high90","high95","mean","sd")
rownames(phi) <- sub("\\phi.", "", rownames(phi))

# Save lambda/psi as one .csv file and save country ideal points as another .csv file
write.csv(lambda,file="lambda_threat9b.csv")
write.csv(phi,file="country_irt_threat9b.csv")



